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TWO  METHODS  FOR  CALCULATING  THE  PHYSICAL  OPTICS  FIELD 


1.  INTRODUCTION 

Given  a  plane  wave  incident  on  a  reflecting  surface  S,  the  physical 
optics  (PO)  field  scattered  by  S  is  given  by 

(1)  J  »  -  -  f  exp[-i  k(  *x]  (  £-u>)  »n  dA 

4  irR  Js 

where  k  »  2  it/  A  is  the  wave  number,  R  is  the  radar  range,  £  is  a  unit  vector 
perpendicular  to  the  wave  front  pointing  from  the  source  to  S,  u)  is  the  unit 
vector  pointing  from  S  to  the  field  point  (observer),  x  is  the  position  vec¬ 
tor  of  a  general  point  on  S,  n  is  the  unit  normal  to  S,  and  dA  is  the  infini¬ 
tesimal  element  of  area  on  S.  Hence  J  is  a  dimensionless  quantity,  and  if  P^ 
is  the  incident  field  power  on  S,  the  power  at  the  field  point  is  given  by  P 

-  pJjI2- 

We  shall  hereafter  be  concerned  with  the  monostatic  case  5  *  -  u>,  so 
that  ignoring  the  range  dependence  in  (1),  the  integral  becomes 

(2)  J  =  ( 1/  A)  1  exp(i  2k  a)  •  x]  to  •  n  dA  . 

JS 

The  radar  cross  section  a  is  then  related  to  J  by 

(3)  a  =■  4ir  |j  |2  . 

The  integral  (2)  is  always  taken  over  only  the  illuminated  part  of  a  surface, 
defined  by  the  relation  u)  •  n  >_  0. 

Closed  form  expressions  for  (2)  can  be  obtained  in  only  a  few  very 
special  cases,  and  the  numerical  evaluation  of  these  integrals  is  difficult. 
However,  in  previous  work  [1]  the  author  has  shown  that  the  double  integral 
(2)  over  S  can  be  reduced  to  a  line  integral  over  the  boundary  of  S  when  S  is 
flat.  If,  moreover,  S  is  a  polygonal  plate  this  line  integral  can  be  further 
reduced  to  a  close-form  expression  involving  no  integrations  at  all.  The 
existence  of  simple  closed-form  formulas  for  polygonal  plates  motivates  the 
facet  decomposition  method  for  calculating  the  PO  field  scattered  by  a  curved 
surface  S:  the  surface  S  is  approximated  by  a  collection  of  triangular 
facets,  and  the  closed  form  scattering  formula  is  applied  to  each  facet.  The 
closed-form  scattering  formulas  for  polygonal  plates  are  reproduced  here  in 
Appendix  A,  with  certain  modifications  which  are  necessary  to  insure  the 
proper  phasing  between  the  various  facets. 

We  shall  also  consider  the  numerical  evaluation  of  (2)  by  means  of  a 
Monte  Carlo  method  in  which  integrals  of  the  form 
Manuscript  approved  January  11,  1984, 


are  approximated  by  sums  of  the  type 

N 

(5)  J  =*  -  ^  f  ( xn) 

N  n=l 

where  the  x^  are  a  random  sample  of  points  from  a  uniform  distribution  on  S. 
In  theory,  there  is  no  limit  on  the  accuracy  of  either  method,  but  there  are 
practical  limits  to  their  usefulness,  and  the  relative  merits  of  the  two 
methods  is  a  question  of  economics  (cost),  and  is  machine  dependent. 

The  Monte  Carlo  method  almost  always  requires  the  use  of  thousands  of 
functional  evaluations,  even  when  S  is  flat,  whereas  the  facet  decomposition 
method  can  sometimes  yield  accurate  results  with  only  a  few  facets  when  S  is 
only  slightly  curved.  On  the  other  hand,  the  Monte  Carlo  method  runs  very 
fast,  the  coding  is  simple,  and  the  requirements  for  computer  storage  are 
minimal,  whereas  the  facet  decomposition  method  requires  a  large  amount  of 
storage.  For  example,  when  the  so  called  "optimal”  compiler  is  used  with  the 
TI  ASC,  one  is  not  permitted  to  use  an  array  with  more  than  2^  -  33,000 
elements,  and  this  limit  is  attained  when  the  number  of  facets  exceeds  5000. 

The  only  practical  limitation  to  the  Monte  Carlo  method  appears  to  be 
round-off  error,  caused  by  finite  word  length.  This  appears  to  limit  the 
number  of  sample  points  to  a  few  tens  of  thousands  for  the  ASC  (which  has  a 
short  word  length),  and  this  limit  is  reached  long  before  one  exceeds  a 
reasonable  cost  limit  for  CPU  time.  The  finite  word  length  also  affects  the 
reliability  of  the  facet  decomposition  method,  but  in  another  way.  The 
scattering  formula  for  triangles  is  of  the  form  (A+B+C)/e,  where,  for 
near  normal  incidence,  e  and  (A+B+C)  are  small,  while  each  of  the  terms  A,  B 
and  C  is  on  the  order  of  unity.  When  this  happens  the  resulting  calculation 
may  be  greatly  in  error,  and  one  has  to  exercise  care  in  devising  methods  to 
detect  the  occurrence  of  such  bad  cases  and  make  the  proper  adjustments.  The 
facet  decomposition  method  can  also  be  expensive  in  CPU  time  since  it 
requires  9  trigonometric  function  evaluations  for  each  facet  whereas  the 
Monte  Carlo  method  requires  only  two  trigonometric  function  evaluations  per 
sample  point. 

Hence,  to  summarize,  the  Monte  Carlo  method  runs  fast,  but  its  accuracy 
is  limited  by  round-off  error.  The  facet  decomposition  method  is  sometimes 
more  costly  in  computer  running  time,  and  is  also  limited  by  constraints  on 
computer  storage.  However,  the  facet  decomposition  method  is  more  accurate 
and  economical  when  the  surface  S  is  only  slightly  curved.  Both  methods 
require  a  substantial  amount  of  calculation  (thousands  of  points  or  facets) 
when  S  is  a  "moderately"  curved  surface  which  contains  only  a  few  Fresnal 
zones. 


Finally,  we  should  remark  that  the  theory  of  Physical  Optics  is  a  high 
frequency  approximation  to  reality,  and  that  the  problem  of  determining  when 


PO  theory  yields  physically  valid  results  is  different  from  the  mathematical 
problem  addressed  in  this  report,  viz.,  the  accurate  numerical  calculation  of 
the  integral  (2). 

2.  TRIANGULAR  DECOMPOSITION 

Letting  x  =  (x^-,  x2 ,  x^)  denote  a  general  position  vector  for  points  on 
S,  the  surface  S  is  defined  by  a  system  of  parametric  equations. 


(6)  x  =  x(u,v) 

where  (u,v)  varies  over  a  domain  D  in  the  Euclidean  plane.  We  usually  take  D 
to  be  the  unit  disk  {u^  +  v^  <.  1},  or  the  unit  square  (0  u  <_  1,  0  ^  v  <_  1}, 
depending  upon  whether  S  has  one  or  two  boundary  curves.  In  the  former  case 
the  boundary  curve  of  S  is  not  necessarily  a  circle  since  (6)  describes  a 
deformation  of  D.  In  the  latter  case  S  is  topologically  equivalent  to  a 
cylinder.  Geometrically,  a  cylinder  is  obtained  by  gluing  together  two 
opposite  sides  of  a  square.  In  analytic  terms,  this  means  that  the  paramet¬ 
ric  equations  for  a  cylinder  satisfy 

x(u,  0)  =  x(u,  1). 

To  obtain  an  approximation  of  S  by  a  collection  of  facets,  we  triangu¬ 
late  D  and  then  lift  the  triangles  via  (6).  Once  a  triangular i zation  is 
given,  a  finer  triangularization  is  obtained  by  replacing  each  triangle  by 
four  smaller  ones,  as  we  now  describe. 

Triangles  in  the  disk  D  are  said  to  be  boundary  or  interior,  the  former 
being  triangles  with  two  vertices  on  the  boundary  circle.  Thus,  in  Figure  1, 
only  triangle  1  is  boundary  (since  2  has  only  one  vertex  on  the  boundary). 

To  subdivide  an  interior  triangle,  we  join  the  midpoints  of  the  legs,  as 
shown  by  the  dotted  lines  in  the  figure.  The  subdivision  of  a  boundary  tri¬ 
angle  is  also  shown  in  Figure  1,  where  A"  is  the  midpoint  of  the  arc  BC. 

Note  that  in  subdividing  a  boundary  triangle,  we  obtain  two  interior  and  two 
boundary  triangles,  and  that  only  one  of  these  four  is  contained  in  the 
original  triangle. 

To  triangulate  the  disk  D  we  first  start  with  an  initial  configuration 
of  F0  triangles  obtained  by  joining  F0  equispaced  points  on  the  boundary 
circle  to  the  center.  Figure  2a  shows  the  initial  configuration  for  the  case 
F0  =  4,  and  it  is  seen  that  all  the  triangles  in  the  initial  configuration 
are  of  boundary  type.  Figure  2b  shows  the  results  of  the  first  subdivision; 
there  are  16  triangles,  and  in  this  case  only  the  odd  numbered  ones  are  of 
the  boundary  type.  At  the  second  subdivision  (not  shown)  there  will  be  64 
triangles,  of  which  only  16  will  be  of  the  boundary  type.  More  generally,  at 
each  subdivision  the  total  number  F  of  triangles  (or  facets)  increases  by  a 
factor  of  4,  while  the  number  of  boundary  triangles  doubles. 

The  initial  configuration  of  a  square  is  always  taken  to  be  the  four 
triangles  shown  in  Figure  3A.  In  this  case  we  do  not  have  to  distinguish 
between  boundary  and  interior  triangles;  at  each  subdivision  each  triangle  is 
decomposed  into  four  smaller  ones  by  joining  the  midpoints. 


V 


Given  a  configuration  of  F  triangles,  let  E,  and  V  denote  the  number 
of  edges  and  vertices,  and  let  Vb  be  the  number  of  vertices  which  lie  on  the 
boundary.  Then,  from  topological  generalities  it  is  known  that 

F  -  E  +  V  =  1, 

2E  =  3F  +  Vb  • 

Hence,  for  F  large,  we  have  the  approximate  relations 

E  »  (3/2)  F  , 

V  «  (1/2)  F  . 

Table  1  shows  the  values  of  F,  E,  V,  Vb  for  up  to  5  subdivisions  and  two 
initial  configurations.  The  row  labled  0  shows  the  data  for  the  initial 
configuration.  Also,  when  D  is  the  disk,  Vb  is  also  the  number  of  boundary 
triangles. 

TABLE  1 

Values  of  F,  E,  V  and  Vb  vs 
Number  of  Iterated  Subdivisions 


Number  of 
Iterations 


Fo  =  3 


F„  =  4 


F 

E 

V 

Vb 

F 

E 

V 

0 

3 

6 

4 

3 

4 

8 

5 

1 

12 

21 

10 

6 

16 

28 

13 

2 

48 

78 

31 

12 

64 

104 

41 

3 

192 

300 

109 

24 

256 

400 

145 

4 

768 

1176 

409 

48 

1024 

1568 

545 

5 

3072 

4656 

1585 

96 

4096 

6208 

2113 

3.  THEORETICAL  ERROR  ANALYSIS 

3A.  ERROR  ANALYSIS  FOR  THE  MONTE  CARLO  METHOD 

To  begin,  given  a  general  integral  of  type  (4)  we  first  calculate  the 
standard  deviation  of  the  estimate  (5).  The  derivation  is  given  in  Appendix 
B,  and  it  turns  out  that  j  is  an  unbiased  estimate  of  J  whose  standard 
deviation  SJ  is  given  by 
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where  A  is  the  area  of  S.  It  is  interesting  to  note  that  from  the  Schwarz 
inequality  we  have 


IJI2  = 


1 


1  dA 


1  /  |£|2  J 


1  dA 


1 


fl2  dA 


with  equality  only  when  f  is  constant;  hence  the  term  in  brackets  measures 
the  variability  of  f. 


In  practice,  the  integral  (2)  ovei  S  is  reduced  to  an  integral  over 
some  planar  domain,  (e.g.,  by  means  of  the  parametric  representation  (6) 
discussed  in  the  previous  section) ,  and  an  especially  revealing  relation  is 
obtained  if  the  plane  P  is  perpendicular  to  the  vector  m.  We  recall  that 
the  scattering  integral  is  performed  over  only  the  illuminated  portion  of  a 
surface  for  which  m  •  q  >  0,  and  we  have 


dAo  ■  ( to  •  q)  dA 

where  dA0  is  the  projection  of  dA  onto  P.  Hence  from  (2), 


(8) 


I 


(1/X)  exp[i2k  m  •  x]  dAc 


where  SQ  is  the  projection  of  S  onto  P.  Applying  (7)  to  the  integral  (8), 
we  have  |f|  *  1/X,  and  therefore 


1  (9) 

1 

&J  =  — 

r  a2  i 

~  -  |j|2 

/N 

L  X2  J 

"  1/2 
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where  Aq  is  the  area  of  SQ.  Setting 

(10)  aQ  -  Air  A2/X2  , 

from  (8)  we  have  |j|  <_  Aq/ X,  and  from  (3)  we  have 

(11)  a  <_  aQ  , 

and  from  (9)  we  have 

Joq-0 


(12) 


&J  = 


/4ir  N 

We  note  that  a  *  c0  only  when  S  is  flat  and  the  incidence  is  normal. 

Following  (3)  it  is  natural  to  estimate  the  radar  cross  section  o  by  o 


where 

(13)  a  =  4ir  |j|2  , 

In  Appendix  B  it  is  shown  that  a  has  a  bias  B  given  by 


rr*;  •r-Lr^w  o.avvvv'V'-'.''  '-.•  -xv w~; 1  r-- ■  .  • -tv-  - 


v-ti 


(14) 


B(  a)  =  (a0-a)/N  . 


s' • 


The  general  expression  for  the  standard  deviation  5a  of  a  is  extremely  messy; 
however,  when  a  <<  a0  and  the  errors  are  small  we  have  the  approximate 
relation 


(15) 


5a/  a  »  /2/N  / a0/  a  . 


Examples  will  be  given  in  Section  3,  and  it  will  be  seen  that  accurate 
estimates  are  harder  to  attain  when  a  is  small. 

B.  Error  Analysis  for  the  Facet  Decomposition  Method 


Let  Rj_,  R?  be  the  principal  radii  of  curvature  on  S  and  let  R0  be  a 
lower  bound  to  R]_,  R2  which  holds  everywhere  on  S.  In  Appendix  C  it  is 
shown  that  when  the  total  area  of  the  collection  of  facets  approximates  the 
area  of  S,  then  the  error  AJ  in  the  estimate  of  J  satisfies 


(16) 


8  IT 

AJ  <  — 


"  ft 


1 

“  • 
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and  for  small  errors  the  corresponding  bound  for  the  relative  error  Aa/ a  is 


(17) 


32  it3/2  a2 

A  a/  a  <  -  •  - - 

“  /3  X2R0/S 


1 

“*  • 
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In  particular,  when  S  is  flat,  Rq  *»  “>,  and  the  right  hand  sides  of  (16)  and 
(17)  are  zero.  Strictly  speaking,  these  relations  are  derived  only  for  the 
case  when  the  facets  are  equilateral  triangles;  however,  we  have  found  that 
these  relations  are  valid  in  the  more  general  case,  except  when  both  the 
actual  and  predicted  relative  errors  are  extremely  small  (<<  10~^).  However, 
the  error  bound  (17)  is  very  often  overly  pessimistic,  and  we  shall  therefore 
discuss  how  the  accuracy  of  the  Facet  Decomposition  method  can  be  judged  by 
observing  how  the  estimates  vary  with  N. 


4.  EXAMPLES 


In  this  section  we  shall  present  the  results  of  numerical  calculations 
for  two  test  cases,  the  sphere  and  the  spherical  cap.  In  both  cases  we 
shall  compare  the  results  of  the  Monte  Carlo  and  Facet  Decomposition  methods 
with  each  other,  and  with  the  performance  predicted  by  the  error  analysis 
discussed  in  the  previous  section. 


4 A.  The  Sphere 


We  first  consider  the  results  for  the  Monte  Carlo  method.  For  a 
sphere  of  radius  a,  the  illuminated  part  is  a  hemisphere  with  area  A  = 
2 it  a^ ,  and  the  planar  projection  SQ  has  area  AQ  =  tt  a^ .  Hence 

Of  =  (4tt/X2)  (tt  a2)2 


6 


When  ka  >  1,  we  have  a  *  it  a^  <<  aQ ,  and  from  (14)  and  (15)  we  have 

4ir2 

(18)  B  =  -  , 

X2  N 

( 19 )  6a/  a  =  ka  /2/N  . 

In  the  cases  considered  below,  we  take  a  =  1  and  X  =  1,  1/2,  1/4,  1/8.  For 
these  cases  we  have  a  =  ir  exactly,  and  in  Tables  2-5  the  columns  labelled  a" 
are  the  estimated  values  of  the  normalized  res  a'  =  a/ tt a2  (->1)  obtained  by 
a  Monte  Carlo  approximation  with  the  indicated  number  N  of  points,  and  B'  and 
6a/ a  are  calculated  according  to 

4  tt2 

B "  =  -  , 

X2  N 


(  2  VT>  TT 
6a/  a  =  -  — 

x/n 

When  6a/ a  >  1,  the  error  is  large  and  the  results  are  unreliable.  Since 
the  true  value  of  a'  is  unity,  the  relative  error  is  |  a'  -  1 1 ,  and  we  note 
that  the  errors  always  fall  within  the  expected  range.  We  also  note  that 
the  errors  would  not  be  significantly  improved  by  making  the  obvious  bias 
correction. 

We  now  consider  the  results  for  the  Facet  Decomposition  Method.  By 
convention,  we  take  the  parametrizing  disc  D  to  lie  in  the  equatorial  plane, 
and  the  entire  sphere  is  triangulated  by  lifting  the  triangularization  of  the 
disc  onto  both  the  northern  and  southern  hemispheres.  As  the  aspect  angle 
is  varied  different  hemispheres  are  illuminated,  and  hence  different  aspect 
angles  correspond  in  effect  to  different  triangularizations  of  a  hemisphere. 

The  results  of  the  calculation  are  shown  in  Tables  6-9  for  four  differ¬ 
ent  aspect  angles  9,  with  9  being  the  colatitude,  and  three  different 
levels  of  triangularization.  The  inital  configuration  (on  the  disc)  was  four 
triangles,  and  we  show  the  results  for  the  3'rd,  4’th,  and  5'th  iterated 
subdivisions,  with  the  number  N  of  facets  being  256,  1024,  and  4096,  respec¬ 
tively.  (cf.  Table  1.)  The  error  bounds  for  Aa/ a  given  by  (17)  are  not 
useful  for  these  cases  since  they  are  all  greater  than  unity,  except  for  the 
case  X  =  1.0,  N  =*  4096  when  the  value  is  .56. 

Confining  our  attention  to  the  column  N  =  4096,  we  see  that  there  are 
some  particularly  bad  results  at  (  X  =  0.5,  3  =  90°);  (X  =*  0.25,  9  =  90°); 

(X  =0.125,  9=0°).  In  fact,  in  the  first  of  these  cases  the  results  are 
better  at  N  3  1024  than  N  =  4096.)  In  Section  1  we  discussed  a  numerical 
problem  which  can  occur  when  the  incidence  is  almost  normal  to  a  facet,  and 
the  computer  was  programmed  to  treat  near  normal  incidence  as  normal  inci¬ 
dence,  and  to  print  a  warning  if  the  contribution  from  a  facet  was  larger 


TABLES  2-5 


Monte  Carlo  Estimates  of  the  Normalized  RCS  o' 
for  the  Unit  Sphere 


Table  2 

Table  4 

X  =  1.0 

X  =  0.25 

N 

o' 

B" 

So/  a 

N 

o'  B' 

So/  a 

1.0523 

.0329 

.26 

1200 

1.8505  .5264 

>  1 

.8903 

.0082 

.13 

4800 

1.7256  .1316 

.51 

.9519 

.064 

.7201  .0329 

.26 

.9686 

.045 

.8137  .0x64 

.18 

Table  3 

Table  5 

X  =  0.5 

i 

X  =  .125 

N 

o' 

B' 

So/  o 

N 

o'  B' 

So/  o 

1200 

.5779 

.1316 

.51 

1200 

4.6054  2.1055 

>  1 

4800 

1.1409 

.0329 

.26 

4800 

3.3170  .5264 

>  1 

19200 

1.0215 

.0082 

.13 

19200 

1.2680  .1316 

.51 

38400 

1.0794 

.0041 

.091 

38400 

1.1914  .0658 

.36 
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TABLES  6-9 


Facet  Decomposition  Estimates  of  the  Normalized  RCS  a' 
for  the  Unit  Sphere 


Table  6 

Table  8 

X  = 

1.0 

X  = 

.25 

N 

0 

256 

1024 

4096 

N 

0 

256 

1024 

4096 

0° 

1.0485 

1.0826 

1.0442 

0° 

.1441 

.7604 

1.1206 

30 

1.0186 

.9184 

.9402 

30 

.3090 

.7550 

1.0431 

60 

.9534 

.9870 

1.0132 

60 

.8500 

.8029 

.9863 

90 

.2198 

.6261 

.9726 

90 

.6208 

.2749 

.5387 

Table  7 

Table  9 

X  = 

0.5 

X  = 

.125 

N 

0 

256 

1024 

4096 

N 

0 

256 

1024 

4096 

0° 

.6975 

1.1027 

1.0957 

0° 

.0300 

.0245 

.7915 

30 

.7164 

.8782 

.9605 

30 

.2246 

1.0458 

1.0154 

60 

.9127 

.8558 

.9326 

60 

.8  220 

1.0939 

.9315 

90 

.0956 

.9905 

.5570 

90 

1.1214 

.7421 

.9890 

than  theory  permits.  Such  warnings  did  not  occur,  and  since  the  worst  (and 
best)  results  did  not  always  occur  at  the  same  aspect  angle,  we  do  not 
believe  that  this  numerical  problem  was  responsible  for  these  bad  results. 
Because  of  the  existence  of  these  bad  cases,  we  cannot  claim  that  the  Facet 
Decomposition  method  is  superior  to  the  Monte  Carlo  method  for  calculating 
the  field  scattered  by  a  sphere. 

4B.  The  Spherical  Cap 


Given  a  sphere  of  radius  a,  the  surface  area  A  and  the  projected  planar 
area  Aq  are  given  by 

A  =  2irah  , 

AQ  =  irr^  =  it  (2ah-h^) 

(See  Figure  4.)  In  the  example  we  shall  keep  A  and  X  fixed,  with  A  =  2ir 
and  X  =  1/8,  and  we  shall  vary  a  in  such  a  manner  that  h  is  in  integral 
multiple  of  X/4.  Since  A  =  2irah  =  2ir,  we  have 

h  =  m  X/4  =  m/32  , 

a  =  1/h  =  32/m 

We  shall  only  present  the  results  for  axial  incidence  (along  the  axis  of 
symmetry) ,  and  hence  m  is  the  number  of  Fresnel  zones  contained  in  each  cap. 
One  would  therefore  expect  that  a  is  small  when  m  is  even,  and  this  turns 
out  to  be  the  case.  Although  the  Fresnel  zones  within  a  given  cap  have 
equal  areas,  a  is  not  exactly  zero  when  m  is  even  because  of  the  "obliquity 
factor"  cos  9  =  m*n  occurring  in  the  integrand  (2). 

The  results  of  the  numerical  calculations  are  shown  in  Table  10  for 
both  the  Monte  Carlo  method  with  40,000  sample  points  and  the  Facet  Decompo¬ 
sition  method  with  4096  facets.  The  values  of  o  are  shown  in  both  absolute 
terms  and  dB,  and  the  standard  deviations  6o  and  error  bounds  A  a  are  calcu¬ 
lated  according  to  (15)  and  (17),  where  a  is  given  the  estimated  value  a 
(unlike  the  case  for  the  sphere,  where  the  "true"  value  of  a  was  known  and 
used).  The  value  of  aQ  is  calculated  according  to  (10),  so  that  aQ  = 

31750  (45.0  dB).  We  also  recall  that  a0  is  the  limiting  value  of  a  as  a 
increases  without  bound,  and  it  was  observed  that  for  large  values  of  a  (not 
shown  in  the  table),  the  estimated  values  of  a  converged  to  this  value. 

We  note  that  for  m  odd  (and  a  large),  the  Monte  Carlo  and  Facet 
Decomposition  estimates  are  within  a  few  tenths  of  a  dB  of  each  other  (except 
when  m  =  9),  and  that  the  two  values  always  differ  by  less  than  6a  . 

For  even  values  of  m  ( a  small),  it  appears  that  the  Facet  Decomposi¬ 
tion  method  gives  better  results.  The  nulls  predicted  by  this  method  are 
deeper  than  those  predicted  by  the  Monte  Carlo  method,  and  one  would  except 
to  see  extremely  small  values  of  a  when  ra  is  small  and  even.  Moreover,  the 
bias  of  the  Monte  Carlo  method  (calculated  according  to  (14))  is  equal  to 
.79,  which  is  larger  than  any  of  the  estimated  values  of  a  for  m  even,  and 
the  values  of  6a  are  also  larger  than  a,  so  that  the  Monte  Carlo  estimates 
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TABLE  10 


Monte  Carlo  and  Facet  Decomposition  Estimates  of  a  for  a 
Spherical  Cap  with  Radius  a  and  Height  h  =  mX/4.  (X  =  0.125) 


Monte  Carlo  Facet  Decomposition 


N 

=  40,000 

N  =  4096 

m 

a 

a 

a  (dB) 

So 

a 

a  (dB) 

Ao 

1 

32.0 

12814.0 

41.1 

110.3 

12764.0 

41.1 

224.0 

3 

10.7 

1375.5 

31.4 

45.7 

1412.7 

31.5 

222.9 

5 

6.4 

525.8 

27.2 

28.6 

499.9 

27.0 

221.7 

7 

4 .6 

236.2 

23.7 

19.3 

250.4 

24.0 

218.3 

9 

3.6 

124.6 

21.0 

14.0 

146.6 

21.7 

213.5 

11 

2.9 

88.4 

19.5 

11.8 

94.6 

19.8 

212.9 

2 

16.0 

.6441 

-1.9 

1.01 

.0001 

-40.0 

0.04 

4 

8.0 

.7228 

-1.4 

1.07 

.0308 

-15.1 

1.39 

6 

5.3 

.4908 

-3.1 

.88 

.1186 

-9.3 

4.12 

8 

4.0 

.1413 

-8.5 

.47 

.2206 

-6.6 

7.45 

10 

3.2 

.4912 

-3.1 

.88 

.3360 

-4.7 

11.50 

12 

2.7 

.5935 

-2.3 

.97 

.4690 

-3.3 

16.10 

are  unreliable.  We  also  note  that  for  m  even  and  increasing  the  Facet 
Decomposition  estimates  also  increase,  which  is  what  one  would  expect  because 
of  Che  variation  of  the  obliquity  factor. 


We  again  emphasize  a  distinction  between  (15)  and  (17).  The  first  is 
an  equation  for  the  standard  deviation  of  an  estimate  involving  random  fac¬ 
tors,  whereas  the  second  is  an  upper  bound  to  the  error  in  an  estimate 
involving  no  random  quantities.  A  large  value  of  5a/ a  calculated  according 
to  (15)  necessarily  implies  that  the  Monte  Carlo  estimate  is  unreliable, 
whereas  a  large  value  for  the  bound  on  Ac/  a  given  by  (17)  carries  no  such 
implication  for  the  Facet  Decomposition  estimate.  However,  when  the  bound 
for  Aa/a  is  large  the  accuracy  of  the  Facet  Decomposition  method  can  be 
judged  by  observing  the  variation  of  o  with  increasing  N.  In  Table  11  we 
have  tabulated  o  vs.  N  for  4  odd  values  and  4  even  values  of  m.  The  esti¬ 
mates  for  m  »  1,  5,  9  appear  to  be  highly  reliable  since  the  values  for  N  = 
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TABLE  11 


Facet 

Decomposition  Values  of  a 

vs.  N  = 

Number  of  : 

Facets 

a  =  32 

a  =  6.4 

a  =  3.6 

a  =  2.9 

a  =  16 

a  =  5.3  a 

=  3.2  a 

=  2. 

N 

m  *  1 

m  =  5 

m  =  9 

m  =  11 

m  =  2 

m  =  6  m 

=  10  m 

=  12 

4 

7323 

206.8 

59.69 

38.29 

1298.5 

139.8 

47.2 

31.3 

16 

11801 

344.4 

26.74 

57.45 

222.7 

223.9 

46.8 

7.0 

64 

12581 

580.4 

180.17 

134.39 

7.9 

9.2 

16.8 

26.1 

256 

12746 

473.8 

164.44 

112.48 

.4472 

.8577 

4.7 

2.4 

1024 

12795 

498.2 

147.11 

98.51 

.0321 

.0921 

.3730 

.2126 

4096 

12764 

499.9 

146.63 

94.63 

.0001 

.1186 

.3360 

.4690 

a  = 

32/m  decreases  as  m 

increases 

,  and  we 

note  that 

acceptable 

results 

are 

obtained  for  the  largest  value  of  a  (flatest  cap)  with  only  64  facets.  For 
the  even  values  of  ra,  we  note  a  large  percentage  difference  between  the 
values  for  N  =  1024  and  N  =  4096,  and  it  therefore  appears  that  the  relative 
errors  are  large.  However,  in  absolute  terms,  one  can  still  feel  confident 
that  a  <  1  in  these  cases,  and  hence  that  by  introducing  small  (quarter- 
wavelength)  changes  in  the  height  of  a  spherical  cap  can  increase  or  decrease 
its  res  by  orders  of  magnitude. 

5.  SUMMARY  AND  CONCLUSIONS 

(i)  The  accuracy  of  the  Monte  Carlo  method  is  limited  by  roundoff  error, 
whereas  the  Facet  Decomposition  method  is  limited  by  computer  storage  and  CPU 
time.  The  Facet  Decomposition  method  is  superior  when  the  reflecting  surface 
S  is  sufficiently  flat. 

(ii)  The  Monte  Carlo  estimate  of  a  has  a  bias  B  and  standard  deviation  6a 
given  by 

(14)  B  =*  (  a0-o)/N  , 

(15)  6  a/ a*  ^2/N  /  a0/  a 

where  N  is  the  number  of  sample  points,  and 
(10)  oQ  -  4itAq/X2  , 

where  Aq  is  the  area  of  the  projection  of  S  onto  a  plane  perpendicular  to 
the  incidence  vector.  An  upper  bound  to  the  error  of  the  Facet  Decomposition 
method  is  given  by 
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where  N  is  the  number  of  facets  and  R0  is  a  lower  bound  to  the  principal 
radii  of  curvature.  (Thus  R0  =  «  when  S  is  flat,  in  which  case  the  Facet 
Decomposition  method  is  exact.)  According  to  these  formulas,  the  calcula¬ 
tion  of  a  should  be  easier  and  more  reliable  when  o,  X,  and  R0  are  "large" 
and  this  prediction  has  been  borne  out  by  experience. 

(iii)  A  large  value  of  <$o/  a  (calculated  according  to  (15)  necessarily 
implies  that  the  Monte  Carlo  estimate  is  unreliable,  whereas  a  large  value 
for  the  bound  on  Aa/ a  given  by  (17)  carries  no  such  implication  for  the 
Facet  Decomposition  estimate.  However,  when  this  bound  is  large,  one  can 
still  gauge  the  accuracy  of  the  Facet  Decomposition  estimate  by  observing 
how  it  varies  with  increasing  N ,  since  accurate  estimates  will  converge  to 
some  value  (the  true  one)  as  N  increases. 
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Fig.  2  —  Subdivision  of  the  disc 
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FIRST  SUBDIVISION 


4  —  Spherical  cap  with  height  =  h 
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APPENDIX  A 


Scattering  Formula  for  a  Flat  Polygonal  Plate 

In  Ref  [1]  we  gave  a  closed  form  expression  for  the  Physical  Optics 
(PO)  field  scattered  by  a  flat  polygonal  plate,  and  for  convenience  the  point 
of  zero  reference  phase  was  chosen  to  lie  somewhere  on  the  polygon.  In  our 
present  work  we  are  concerned  with  summing  the  contributions  from  many  dif¬ 
ferent  facets,  and  their  relative  phases  must  now  be  taken  into  account.  A 
convenient  way  of  doing  this  is  to  take  the  point  of  zero  reference  phase  to 
be  the  field  point,  and  to  refer  the  vertices  of  all  the  polygons  to  a  common 
coordinate  system.  (It  might  be  "natural"  to  assume  that  the  phase  of  the 
field  scattered  by  each  facet  corresponds  to  its  centroid,  but  this  is  not 
correct.)  To  obtain  the  correct  result,  we  multiply  the  expression  for  the 
incident  field  (given  at  the  bottom  of  p.  590  of  Ref  [1])  by  a  certain  quan¬ 
tity  (exp[-  ik  £  •  R] )  which  makes  the  phase  zero  at  the  field  point.  Then, 
by  linearity,  we  multiply  all  the  expressions  given  in  Ref.  [1]  for  the 
scattered  field  by  the  same  quantity. 

We  shall  give  the  resulting  equations  for  a  flat  polygonal  plate  with 
N  vertices.  The  results  will  be  stated  for  the  bistatic  case,  and  we  use  the 
same  notation  as  in  equation  (1)  of  Section  1  of  this  report.  In  addition, 
let 


6  =  fixed  origin  (same  for  all  facets) 
f  =  position  vector  of  the  field  point  P; 
i.e. ,  f  =  P  -  0. 

xn  =  position  vector  of  the  n'th  vertex  of  the  polygon, 
a  =  £  -  ai  (=  -2oj  for  the  raonostatic  case) 

6  =  a  -  ( a»n)  n  . 


We  set  xN+i  =  xi  and  A  xn  =  xn+i  -  xn,  (n  =  1, 
by  the  polygon  is  given  by 

For  B  *  0, 


(4*  R)  J  = 


exp(-  ika»f) 


where 


sin(y  a*A  xn) 


2,. 


N). 


n=l 


Tn 


The  field  J  scattered 


) 


\  ik  ( 

ex?  j  ~  a’(xn+xn+l)  j  • 


a 


(4ir  R) 


[-  exp(-  ika*f)]  [2ik  (m*n)  A  exp(ika*q)] 


where  A  is  the  area  of  the  polygon  and  q  is  the  position  vector  of  any  point 
on  the  polygon.  (When  3*0,  o  is  a  scalar  multiple  of  the  normal  vector  n, 
and  hence  a*(q-q')  =  0  for  any  two  points  q,  q '  on  the  polygon.) 

We  note  that  the  factor  [-  exp(-  ik  a*f)]  is  the  same  for  all  facets, 
and  hence  can  be  neglected  when  the  total  set  of  scatterers  consists  of 
polygons  only. 
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APPENDIX  B 


Error  Analysis  for  the  Monte  Carlo  Method 

Let  x  be  a  random  point  with  a  uniform  distribution  on  S,  and  let  f  be 
a  function  on  S  with  integral 


J 


Is 


f  dA  . 


The  expected  value  of  f(x)  is 


E  f(x) 


I 


dA 

f(x)  —  =  j/A 
A 


since  dA/A  is  the  uniform  density  on  S.  Hence  the  quantity  A  f(x)  is  an 
unbiased  estimator  of  J,  and  its  variance  is 


E[  |  A  f(x)|2] 


|J| 


=  J  A|f|2  dA  -  | J | 2  . 


The  relation  (7)  then  follows  from  (5),  since  the  quantities  A  f(xn)  are 
uncorrelated  estimates  of  J. 

The  calculation  of  the  standard  deviation  &J  is  given  in  the  text. 

To  calculate  the  bias  B,  let  AJ  be  the  error  in  j.  Then 

a  =  4-rr  I J  1 2  «  4  ir  |j  +  Aj|2  . 

Expanding,  and  using  the  fact  that  E(AJ)  =  0  (since  J  is  unbiased),  we  get 
Eo=  4tt[  ( J  |  2  +  E  |  Aj|2]  =  a  +  4tt  |  6j  j 2  . 


Hence,  using  (12), 


E  a  -  a  =  4  tt  |  6j  |  2  =  (  a  -  a)  /  N  . 


We  conclude  by  sketching  the  derivation  of  (15).  Let  and  J2  be  the 
real  imaginary  parts  of  J,  so  that 


4ir(J12  +  J22) 


Then,  neglecting  higher  order  terms  for  small  errors,  we  have 

Ao  =  8 it  (j]_  AJ]_  +  J2  AJ2)  • 


•_>  'j  • .  •••■ 


;•  .• 
w*  ■_ 
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Squaring,  and  applying  che  expected  value  operator  E,  we  get 


APPENDIX  C 


Derivation  of  the  Error  Bounds  for  Facet  Decomposition 
Let  f  denote  the  vector  valued  function 
f(x)  =  exp(i2ku)*x]  a)/  \  , 


so  that 


Since  ju)|2  *  1,  the  divergence  of  f  is  given  by 
V*f  =  ( i2 k./  A)  exp[i2km*x]  . 
Hence,  using  the  Divergence  Theorem, 


All 


7*f  dV 


<  (4tt/A2)  AV  , 


where  AV  is  the  volume  contained  between  the  surface  S  and  the  approximating 
collection  of  facets. 


With  a  little  geometry  one  can  show  that 
AV  <  h0  A  , 

where  hQ  is  an  upperbound  to  the  height  of  any  point  on  S  above  the  nearest 
facet,  and  that 


ho  1 


D2 

2Rr 


where  D  is  an  upperbound  to  the  length  of  the  facet  edges.  Hence 

2itA 


AJ  < 


A2R, 


D2 


If  we  now  assume  that  the  tr iangularization  is  so  fine  that  the  total  area 
of  the  facets  is  approximately  A,  and  that  the  triangles  are  equilateral, 
we  have 


N(#  °2) 


From  these  last  two  relations  we  get  (16). 

To  derive  (17)  from  (16)  we  use  the  usual  first  order  approximations 


